/****************************************************************************
*                                                                           *
* scatmatrix.cpp								                            *
*                                                                           *
* ALaRI 2008												                *
*                                                                           *
****************************************************************************/

/****************************************************************************
* Includes										                            *
****************************************************************************/
#include "scatmatrix.h" // - me
#include <memory>

/****************************************************************************
* CScatMatrixL0 - private class for scattering matrix (level 0)				*
****************************************************************************/

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatMatrixL0::getWidth() const
{
	// only one pixel in simulation space
	return 1;
}

int CScatMatrixL0::getHeight() const
{
	// only one pixel in simulation space
	return 1;
}

const CScatBlock* CScatMatrixL0::getBlock(eBlock nBlock, bool /*bCreateIfNotThere*/) const
{

	// switch the block needed
	switch (nBlock)
	{
	// diagonal
	case eBlock_N2N:
	case eBlock_E2E:
	case eBlock_S2S:
	case eBlock_W2W:
	return &m_fltThetta02;
	// non-dialogal
	case eBlock_N2E:
	case eBlock_N2S:
	case eBlock_N2W:
	case eBlock_E2N:
	case eBlock_E2S:
	case eBlock_E2W:
	case eBlock_S2N:
	case eBlock_S2E:
	case eBlock_S2W:
	case eBlock_W2N:
	case eBlock_W2E:
	case eBlock_W2S:
	return &m_fltThetta01;
	// these matrices don't have I matrix
	case eBlock_ILT:
	case eBlock_IRT:
	case eBlock_ILB:
	case eBlock_IRB:
	break;
	default: // invalid
	assert(FALSE);
	}
	// no block found
	return NULL;
}

const CScatMatrix* CScatMatrixL0::getLOrT() const
{
	// no child matrices
	return NULL;
}

const CScatMatrix* CScatMatrixL0::getROrB() const
{
	// no child matrices
	return NULL;
}

const CScatMatrix* CScatMatrixL0::getTransformed(
	CTransform::eTransform /*nTransform*/, bool /*bCreateIfNotThere*/) const
{
	// simple matrix does not need rotations/reflections in space
	return this;
}

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatMatrixL0::CScatMatrixL0(const std::complex<float>& fltThetta01,
		const std::complex<float>& fltThetta02) :
	m_fltThetta01(fltThetta01),
	m_fltThetta02(fltThetta02)
{
	// nothing to do
}

/****************************************************************************
* CScatMatrixLP - scattering matrix (level+)								*
****************************************************************************/

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatMatrixLP::getWidth() const
{
	return m_nWidth;
}

int CScatMatrixLP::getHeight() const
{
	return m_nHeight;
}

const CScatMatrix* CScatMatrixLP::getLOrT() const
{
	return m_pLOrT;
}

const CScatMatrix* CScatMatrixLP::getROrB() const
{
	return m_pROrB;
}

const CScatBlock* CScatMatrixLP::getBlock(eBlock nBlock, bool bCreateIfNotThere) const
{
	// no substitution of the block assumed by default
	CScatMatrix::eBlock nBlock2Return = nBlock;

	// look at the block requested
	switch (nBlock)
	{
	case eBlock_N2N:
	case eBlock_N2E:
	case eBlock_N2S:
	case eBlock_N2W:
	case eBlock_E2E:
	case eBlock_E2S:
	case eBlock_E2W:
	case eBlock_S2S:
	case eBlock_S2W:
	case eBlock_W2W:
	case eBlock_ILT:
	case eBlock_IRB:
	case eBlock_IRT:
	break;
	// all others are substituted
	case eBlock_E2N:
		nBlock2Return = CScatMatrix::eBlock_N2E;
	break;
	case eBlock_S2N:
		nBlock2Return = CScatMatrix::eBlock_N2S;
	break;
	case eBlock_S2E:
		nBlock2Return = CScatMatrix::eBlock_E2S;
	break;
	case eBlock_W2N:
		nBlock2Return = CScatMatrix::eBlock_N2W;
	break;
	case eBlock_W2E:
		nBlock2Return = CScatMatrix::eBlock_E2W;
	break;
	case eBlock_W2S:
		nBlock2Return = CScatMatrix::eBlock_S2W;
	break;
	case eBlock_ILB:
		nBlock2Return = CScatMatrix::eBlock_IRT;
	break;
	default: // invalid
		assert(FALSE);
	return NULL;
	}

	// first construct the block
	if (m_prgnBlockLPs[nBlock2Return] == NULL && bCreateIfNotThere)
	{
		// construct blocks on demand since not all of them may even be needed
		((CScatMatrixLP*)this)->constructBlockLP(nBlock2Return);
		assert(m_prgnBlockLPs[nBlock2Return] != NULL); // should be there
	}

	// if we substituted we need to transpose - fr270
	if (m_prgnBlockLPs[nBlock2Return] != NULL && nBlock2Return != nBlock)
		return m_prgnBlockLPs[nBlock2Return]->getTransformed(
			CTransform::eTransform_fr270, bCreateIfNotThere);

	// result
	return m_prgnBlockLPs[nBlock2Return];
}

const CScatMatrix* CScatMatrixLP::getTransformed(
	CTransform::eTransform nTransform, bool bCreateIfNotThere) const
{
	assert(nTransform >= CTransform::eTransform_r0 && nTransform < CTransform::eTransform_count);

	// check null transformation first
	if (nTransform == CTransform::eTransform_r0)
		return this;
	// create if not created
	else if (m_rgnpTransforms[nTransform - 1] == NULL && bCreateIfNotThere)
	{
		CScatMatrixLP* pThis = (CScatMatrixLP*)this;
		pThis->m_rgnpTransforms[nTransform - 1] = new CScatMatrixLPT(pThis, nTransform);
	}

	return m_rgnpTransforms[nTransform - 1];
}

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatMatrixLP::CScatMatrixLP(const CScatMatrix *pLOrT,
		const CScatMatrix *pROrB, bool bHorizontal) :
	m_pLOrT(pLOrT),
	m_pROrB(pROrB)
{
	assert(m_pLOrT->getWidth() >= 1 && m_pLOrT->getHeight() >= 1);
	assert(m_pROrB->getWidth() >= 1 && m_pROrB->getHeight() >= 1);
	assert(bHorizontal ?
		(pLOrT->getHeight() == pROrB->getHeight()) : (pLOrT->getWidth() == pROrB->getWidth()));

	// get width and height
	if (bHorizontal)
	{
		m_nWidth = pLOrT->getWidth() + pROrB->getWidth();
		m_nHeight = pLOrT->getHeight();
	}
	else
	{
		m_nWidth = pLOrT->getWidth();
		m_nHeight = pLOrT->getHeight() + pROrB->getHeight();
	}

	// cleanup other data
	memset(&m_rgnpTransforms[0], 0, sizeof(m_rgnpTransforms));
	memset(&m_prgnBlockLPs[0], 0, sizeof(m_prgnBlockLPs));
}

CScatMatrixLP::~CScatMatrixLP()
{
	// delete transforms
	for (int nTransform = CountOf(m_rgnpTransforms); --nTransform >= 0; )
		if (m_rgnpTransforms[nTransform] != NULL)
			delete m_rgnpTransforms[nTransform];

	// delete blocks
	for (int nBlock = CountOf(m_prgnBlockLPs); --nBlock >= 0; )
		if (m_prgnBlockLPs[nBlock] != NULL)
			delete m_prgnBlockLPs[nBlock];
}

/****************************************************************************
* Protected helpers															*
****************************************************************************/
void CScatMatrixLP::constructBlockLP(CScatMatrix::eBlock nBlock)
{
	switch (nBlock)
	{
	// exchange scattering
	// NOTE: we compute several blocks at once because it is faster!
	case eBlock_N2N:
	case eBlock_N2E:
	case eBlock_N2S:
	case eBlock_N2W:
		computeBlockN2X();
	break;
	case eBlock_E2E:
	case eBlock_E2S:
	case eBlock_E2W:
		computeBlockE2X();
	break;
	case eBlock_S2S:
	case eBlock_S2W:
		computeBlockS2X();
	break;
	case eBlock_W2W:
		computeBlockW2W();
	break;
	// I matrix
	case eBlock_ILT:
		computeBlockILT();
	break;
	case eBlock_IRB:
		computeBlockIRB();
	break;
	case eBlock_IRT:
		computeBlockIRT();
	break;
	default:
		assert(FALSE);
	}
}

void CScatMatrixLP::computeBlockN2X()
{
	// test if we are horizontal
	if (getHeight() == m_pLOrT->getHeight())
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spILT_INE(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_N2E, true)));
		std::auto_ptr<CScatBlockLP> spILB_INE(simpleMultiply(
			getBlock(eBlock_ILB, true), m_pLOrT->getBlock(eBlock_N2E, true)));
		std::auto_ptr<CScatBlockLP> spIRT_JNW(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_N2W, true)));
		std::auto_ptr<CScatBlockLP> spIRB_JNW(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_N2W, true)));

		// call helpers to create blocks
		m_prgnBlockLPs[eBlock_N2N] = computeSymetricBlock2x2(m_pLOrT->getBlock(eBlock_N2N, true),
			m_pLOrT->getBlock(eBlock_E2N, true), spILT_INE.get(), spIRT_JNW.get(),
			m_pROrB->getBlock(eBlock_N2N, true), m_pROrB->getBlock(eBlock_W2N, true), spIRB_JNW.get());

		m_prgnBlockLPs[eBlock_N2E] = computeNormalBlock2x1(
			m_pROrB->getBlock(eBlock_W2E, true), spILB_INE.get(),
			m_pROrB->getBlock(eBlock_N2E, true), spIRB_JNW.get());
			
		m_prgnBlockLPs[eBlock_N2S] = computeNormalBlock2x2(
			m_pLOrT->getBlock(eBlock_N2S, true), m_pLOrT->getBlock(eBlock_E2S, true),
			spILT_INE.get(), spIRT_JNW.get(), m_pROrB->getBlock(eBlock_W2S, true), spILB_INE.get(),
			m_pROrB->getBlock(eBlock_N2S, true), spIRB_JNW.get());

		m_prgnBlockLPs[eBlock_N2W] = computeNormalBlock2x1F(m_pLOrT->getBlock(eBlock_N2W, true),
			m_pLOrT->getBlock(eBlock_E2W, true), spILT_INE.get(), spIRT_JNW.get());
	}
	else
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spILT_INS(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_N2S, true)));
		std::auto_ptr<CScatBlockLP> spILB_INS(simpleMultiply(
			getBlock(eBlock_ILB, true), m_pLOrT->getBlock(eBlock_N2S, true)));

		// now call helpers to create blocks
		m_prgnBlockLPs[eBlock_N2N] = computeSymetricBlock1x1(
			m_pLOrT->getBlock(eBlock_N2N, true), m_pLOrT->getBlock(eBlock_S2N, true), spILT_INS.get());
	
		m_prgnBlockLPs[eBlock_N2E] = computeNormalBlock1x2(
			m_pLOrT->getBlock(eBlock_N2E, true), m_pLOrT->getBlock(eBlock_S2E, true),
			spILT_INS.get(), m_pROrB->getBlock(eBlock_N2E, true), spILB_INS.get());		

		m_prgnBlockLPs[eBlock_N2S] = simpleMultiply(
			m_pROrB->getBlock(eBlock_N2S, true), spILB_INS.get());

		m_prgnBlockLPs[eBlock_N2W] = computeNormalBlock1x2(
			m_pLOrT->getBlock(eBlock_N2W, true), m_pLOrT->getBlock(eBlock_S2W, true),
			spILT_INS.get(), m_pROrB->getBlock(eBlock_N2W, true), spILB_INS.get());
	}
}

void CScatMatrixLP::computeBlockE2X()
{
	// test if we are horizontal
	if (getHeight() == m_pLOrT->getHeight())
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spIRT_JEW(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_E2W, true)));
		std::auto_ptr<CScatBlockLP> spIRB_JEW(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_E2W, true)));

		// call helpers to create blocks
		m_prgnBlockLPs[eBlock_E2E] = computeSymetricBlock1x1(
			m_pROrB->getBlock(eBlock_E2E, true), m_pROrB->getBlock(eBlock_W2E, true), spIRB_JEW.get());

		m_prgnBlockLPs[eBlock_E2S] = computeNormalBlock1x2F(
			m_pLOrT->getBlock(eBlock_E2S, true), spIRT_JEW.get(), m_pROrB->getBlock(eBlock_E2S, true),
			m_pROrB->getBlock(eBlock_W2S, true), spIRB_JEW.get());

		m_prgnBlockLPs[eBlock_E2W] = simpleMultiply(
			m_pLOrT->getBlock(eBlock_E2W, true), spIRT_JEW.get());
	}
	else
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spILT_IES(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_E2S, true)));
		std::auto_ptr<CScatBlockLP> spILB_IES(simpleMultiply(
			getBlock(eBlock_ILB, true), m_pLOrT->getBlock(eBlock_E2S, true)));
		std::auto_ptr<CScatBlockLP> spIRT_JEN(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_E2N, true)));
		std::auto_ptr<CScatBlockLP> spIRB_JEN(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_E2N, true)));

		// now call helpers to create blocks
		m_prgnBlockLPs[eBlock_E2E] = computeSymetricBlock2x2(
			m_pLOrT->getBlock(eBlock_E2E, true), m_pLOrT->getBlock(eBlock_S2E, true),
			spILT_IES.get(), spIRT_JEN.get(), m_pROrB->getBlock(eBlock_E2E, true),
			m_pROrB->getBlock(eBlock_N2E, true), spIRB_JEN.get());

		m_prgnBlockLPs[eBlock_E2S] = computeNormalBlock2x1(
			m_pROrB->getBlock(eBlock_N2S, true), spILB_IES.get(),
			m_pROrB->getBlock(eBlock_E2S, true), spIRB_JEN.get());

		m_prgnBlockLPs[eBlock_E2W] = computeNormalBlock2x2(
			m_pLOrT->getBlock(eBlock_E2W, true), m_pLOrT->getBlock(eBlock_S2W, true),
			spILT_IES.get(), spIRT_JEN.get(), m_pROrB->getBlock(eBlock_N2W, true), spILB_IES.get(),
			m_pROrB->getBlock(eBlock_E2W, true), spIRB_JEN.get());
	}
}

void CScatMatrixLP::computeBlockS2X()
{
	// test if we are horizontal
	if (getHeight() == m_pLOrT->getHeight())
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spILT_ISE(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_S2E, true)));
		std::auto_ptr<CScatBlockLP> spILB_ISE(simpleMultiply(
			getBlock(eBlock_ILB, true), m_pLOrT->getBlock(eBlock_S2E, true)));
		std::auto_ptr<CScatBlockLP> spIRT_JSW(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_S2W, true)));
		std::auto_ptr<CScatBlockLP> spIRB_JSW(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_S2W, true)));

		// call helpers to create blocks
		m_prgnBlockLPs[eBlock_S2S] = computeSymetricBlock2x2(m_pLOrT->getBlock(eBlock_S2S, true),
			m_pLOrT->getBlock(eBlock_E2S, true), spILT_ISE.get(), spIRT_JSW.get(),
			m_pROrB->getBlock(eBlock_S2S, true), m_pROrB->getBlock(eBlock_W2S, true), spIRB_JSW.get());

		m_prgnBlockLPs[eBlock_S2W] = computeNormalBlock2x1F(m_pLOrT->getBlock(eBlock_S2W, true),
			m_pLOrT->getBlock(eBlock_E2W, true), spILT_ISE.get(), spIRT_JSW.get());
	}
	else
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spIRB_JSN(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_S2N, true)));
		std::auto_ptr<CScatBlockLP> spIRT_JSN(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_S2N, true)));

		// now call helpers to create blocks
		m_prgnBlockLPs[eBlock_S2S] = computeSymetricBlock1x1(
			m_pROrB->getBlock(eBlock_S2S, true), m_pROrB->getBlock(eBlock_N2S, true), spIRB_JSN.get());

		m_prgnBlockLPs[eBlock_S2W] = computeNormalBlock1x2F(
			m_pLOrT->getBlock(eBlock_S2W, true), spIRT_JSN.get(), m_pROrB->getBlock(eBlock_S2W, true),
			m_pROrB->getBlock(eBlock_N2W, true), spIRB_JSN.get());
	}
}

void CScatMatrixLP::computeBlockW2W()
{
	// test if we are horizontal
	if (getHeight() == m_pLOrT->getHeight())
	{
		// create common block
		std::auto_ptr<CScatBlockLP> spILT_IWE(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_W2E, true)));

		// call helpers to create blocks
		m_prgnBlockLPs[eBlock_W2W] = computeSymetricBlock1x1(
			m_pLOrT->getBlock(eBlock_W2W, true), m_pLOrT->getBlock(eBlock_E2W, true), spILT_IWE.get());
	}
	else
	{
		// create common blocks
		std::auto_ptr<CScatBlockLP> spILT_IWS(simpleMultiply(
			getBlock(eBlock_ILT, true), m_pLOrT->getBlock(eBlock_W2S, true)));
		std::auto_ptr<CScatBlockLP> spIRT_JWN(simpleMultiply(
			getBlock(eBlock_IRT, true), m_pROrB->getBlock(eBlock_W2N, true)));
		std::auto_ptr<CScatBlockLP> spIRB_JWN(simpleMultiply(
			getBlock(eBlock_IRB, true), m_pROrB->getBlock(eBlock_W2N, true)));

		// now call helpers to create blocks
		m_prgnBlockLPs[eBlock_W2W] = computeSymetricBlock2x2(
			m_pLOrT->getBlock(eBlock_W2W, true), m_pLOrT->getBlock(eBlock_S2W, true),
			spILT_IWS.get(), spIRT_JWN.get(), m_pROrB->getBlock(eBlock_W2W, true),
			m_pROrB->getBlock(eBlock_N2W, true), spIRB_JWN.get());
	}
}

void CScatMatrixLP::computeBlockILT()
{
	assert(m_prgnBlockLPs[eBlock_ILT] == NULL);

	const CScatBlock *pBlockROrB = m_pROrB->getBlock(
		getHeight() == m_pLOrT->getHeight() ? eBlock_W2W : eBlock_N2N, true);
	const CScatBlock *pBlockIRT = getBlock(eBlock_IRT, true);

	m_prgnBlockLPs[eBlock_ILT] = symetricMultiply(pBlockIRT, pBlockROrB);
}

void CScatMatrixLP::computeBlockIRB()
{
	assert(m_prgnBlockLPs[eBlock_IRB] == NULL);

	const CScatBlock *pBlockLOrT = m_pLOrT->getBlock(
		getHeight() == m_pLOrT->getHeight() ? eBlock_E2E : eBlock_S2S, true);
	const CScatBlock *pBlockILB = getBlock(eBlock_ILB, true);

	m_prgnBlockLPs[eBlock_IRB] = symetricMultiply(pBlockILB, pBlockLOrT);
}

void CScatMatrixLP::computeBlockIRT()
{
	assert(m_prgnBlockLPs[eBlock_IRT] == NULL);

	// blocks we need from children
	const CScatBlock *pBlockLOrT, *pBlockROrB;

	// test if we are horizontal
	if (getHeight() == m_pLOrT->getHeight())
		pBlockLOrT = m_pLOrT->getBlock(eBlock_E2E, true),
		pBlockROrB = m_pROrB->getBlock(eBlock_W2W, true);
	else
		pBlockLOrT = m_pLOrT->getBlock(eBlock_S2S, true),
		pBlockROrB = m_pROrB->getBlock(eBlock_N2N, true);

	// blocks should be square and of equal size
	int nSize = pBlockROrB->getWidth();
	assert(pBlockROrB->getHeight() == nSize);
	assert(pBlockLOrT->getWidth() == nSize && pBlockLOrT->getHeight() == nSize);

	// now create a new block of apropriate size
	CScatBlockLP *pBlock = new CScatBlockLP(false, nSize, nSize);

	int nRow, nColumn, nRover; // now multiply the two blocks (-sign)
	for (nRow = 0; nRow < nSize; ++nRow)
		for (nColumn = 0; nColumn < nSize; ++nColumn)
			for (nRover = 0; nRover < nSize; ++nRover)
				pBlock->valueAt(nRow, nColumn) -=
					pBlockROrB->valueAt(nRow, nRover) * pBlockLOrT->valueAt(nRover, nColumn);

	// add I to the matrix
	for (nRover = 0; nRover < nSize; ++nRover)
		pBlock->valueAt(nRover, nRover) += 1.f;

	// and invert
	try
	{
		pBlock->invert();
	}
	catch(...)
	{
		// cleanup
		delete pBlock;
	}

	// set new block
	m_prgnBlockLPs[eBlock_IRT] = pBlock;
}

CScatBlockLP* CScatMatrixLP::symetricMultiply(const CScatBlock *pL, const CScatBlock *pR)
{
	// blocks should be square and of equal size
	int nSize = pL->getWidth();
	assert(pL->getHeight() == nSize);
	assert(pR->getWidth() == nSize && pR->getHeight() == nSize);

	// now create a new block of apropriate size
	CScatBlockLP *pBlock = new CScatBlockLP(true, nSize, nSize);

	int nRow, nColumn, nRover; // now multiply two blocks
	for (nRow = 0; nRow < nSize; ++nRow)
		for (nColumn = nRow; nColumn < nSize; ++nColumn)
			for (nRover = 0; nRover < nSize; ++nRover)
				pBlock->valueAt(nRow, nColumn) +=
					pL->valueAt(nRow, nRover) * pR->valueAt(nRover, nColumn);

	// result
	return pBlock;
}

CScatBlockLP* CScatMatrixLP::simpleMultiply(const CScatBlock *pL, const CScatBlock *pR)
{
	// blocks should be of appropriate size
	assert(pL->getWidth() == pR->getHeight());

	// now create a new block of result size
	CScatBlockLP *pBlock = new CScatBlockLP(false, pR->getWidth(), pL->getHeight());

	int nRow, nColumn, nRover; // now multiply the two blocks
	for (nRow = pBlock->getHeight(); --nRow >= 0; )
		for (nColumn = pBlock->getWidth(); --nColumn >= 0; )
			for (nRover = pL->getWidth(); --nRover >= 0; )
				pBlock->valueAt(nRow, nColumn) +=
					pL->valueAt(nRow, nRover) * pR->valueAt(nRover, nColumn);

	return pBlock;
}

CScatBlockLP* CScatMatrixLP::computeSymetricBlock2x2(const CScatBlock *pB1,
	const CScatBlock *pB2, const CScatBlock *pB3, const CScatBlock *pB4,
	const CScatBlock *pB5, const CScatBlock *pB6, const CScatBlock *pB7)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL);
	assert(pB5 != NULL && pB6 != NULL && pB7 != NULL);
	// first expression
	assert(pB1->getWidth() == pB3->getWidth());
	assert(pB2->getWidth() == pB3->getHeight());
	assert(pB1->getHeight() == pB2->getHeight());
	// second expression
	assert(pB2->getWidth() == pB4->getHeight());
	// third expression
	assert(pB5->getWidth() == pB7->getWidth());
	assert(pB6->getWidth() == pB7->getHeight());
	assert(pB5->getHeight() == pB6->getHeight());
	// all three expressions
	assert(pB5->getWidth() == pB4->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(true,
		pB1->getWidth() + pB5->getWidth(), pB1->getHeight() + pB5->getHeight());

	// multiply blocks 2 and 4
	int nRow, nColumn, nRover;
	for (nRow = pB2->getHeight(); --nRow >= 0; )
		for (nColumn = pB4->getWidth(); --nColumn >= 0; )
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, pB1->getWidth() + nColumn) +=
					pB2->valueAt(nRow, nRover) * pB4->valueAt(nRover, nColumn);

	// now what we need is to multiply B2 and B3 and add block 1.
	// but note that we need to do it only for the elements above or on the diagonal.
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = nRow; nColumn < pB1->getWidth(); ++nColumn)
		{
			pB->valueAt(nRow, nColumn) = pB1->valueAt(nRow, nColumn);
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB2->valueAt(nRow, nRover) * pB3->valueAt(nRover, nColumn);
		}

	// B6 and B7 with B5
	for (nRow = pB5->getHeight(); --nRow >= 0; )
		for (nColumn = nRow; nColumn < pB5->getWidth(); ++nColumn)
		{
			pB->valueAt(pB1->getHeight() + nRow, pB1->getWidth() + nColumn) =
				pB5->valueAt(nRow, nColumn);
			for (nRover = pB6->getWidth(); --nRover >= 0; )
				pB->valueAt(pB1->getHeight() + nRow, pB1->getWidth() + nColumn) +=
					pB6->valueAt(nRow, nRover) * pB7->valueAt(nRover, nColumn);
		}

	// result
	return pB;
}

CScatBlockLP* CScatMatrixLP::computeSymetricBlock1x1(
	const CScatBlock *pB1, const CScatBlock *pB2, const CScatBlock *pB3)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL);
	// expression is valid
	assert(pB2->getWidth() == pB3->getHeight());
	assert(pB1->getHeight() == pB2->getHeight());
	assert(pB1->getWidth() == pB3->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(true, pB1->getWidth(), pB1->getHeight());

	// multipy blocks 2 and 3 and add B1 writing results inplace
	int nRow, nColumn, nRover;
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = nRow; nColumn < pB1->getWidth(); ++nColumn)
		{
			pB->valueAt(nRow, nColumn) = pB1->valueAt(nRow, nColumn);
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB2->valueAt(nRow, nRover) * pB3->valueAt(nRover, nColumn);
		}

	// result
	return pB;
}

CScatBlockLP* CScatMatrixLP::computeNormalBlock2x1(const CScatBlock *pB1,
	const CScatBlock *pB2, const CScatBlock *pB3, const CScatBlock *pB4)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL);
	// first expression is valid
	assert(pB1->getWidth() == pB2->getHeight());
	// second expression is valid
	assert(pB1->getWidth() == pB4->getHeight());
	assert(pB3->getHeight() == pB1->getHeight());
	assert(pB3->getWidth() == pB4->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(false,
		pB2->getWidth() + pB3->getWidth(), pB3->getHeight());

	// multiply B1, B4 and add B3 writing results inplace
	int nRow, nColumn, nRover;
	for (nRow = pB3->getHeight(); --nRow >= 0; )
		for (nColumn = pB3->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow, nColumn + pB2->getWidth()) = pB3->valueAt(nRow, nColumn);
			for (nRover = pB1->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn + pB2->getWidth()) +=
					pB1->valueAt(nRow, nRover) * pB4->valueAt(nRover, nColumn);
		}

	// multiply B1, B2 writing results inplace
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = pB2->getWidth(); --nColumn >= 0; )
			for (nRover = pB1->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB1->valueAt(nRow, nRover) * pB2->valueAt(nRover, nColumn);

	// result
	return pB;
}

CScatBlockLP* CScatMatrixLP::computeNormalBlock2x1F(const CScatBlock *pB1,
	const CScatBlock *pB2, const CScatBlock *pB3, const CScatBlock *pB4)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL);
	// first expression is valid
	assert(pB1->getWidth() == pB3->getWidth());
	assert(pB1->getHeight() == pB2->getHeight());
	assert(pB2->getWidth() == pB3->getHeight());
	// second expression is valid
	assert(pB2->getWidth() == pB4->getHeight());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(false,
		pB1->getWidth() + pB4->getWidth(), pB1->getHeight());

	// multiply B2, B3 and add B1 writing results inplace
	int nRow, nColumn, nRover;
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = pB1->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow, nColumn) = pB1->valueAt(nRow, nColumn);
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB2->valueAt(nRow, nRover) * pB3->valueAt(nRover, nColumn);
		}

	// multiply B2, B4 writing results inplace
	for (nRow = pB2->getHeight(); --nRow >= 0; )
		for (nColumn = pB4->getWidth(); --nColumn >= 0; )
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn + pB1->getWidth()) +=
					pB2->valueAt(nRow, nRover) * pB4->valueAt(nRover, nColumn);

	return pB;
}

CScatBlockLP* CScatMatrixLP::computeNormalBlock1x2(const CScatBlock *pB1,
	const CScatBlock *pB2, const CScatBlock *pB3, const CScatBlock *pB4, const CScatBlock *pB5)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL && pB5 != NULL);
	// first expression is valid
	assert(pB1->getWidth() == pB3->getWidth());
	assert(pB1->getHeight() == pB2->getHeight());
	assert(pB2->getWidth() == pB3->getHeight());
	// second expression is valid
	assert(pB4->getWidth() == pB5->getHeight());
	// both expressions
	assert(pB1->getWidth() == pB5->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(false,
		pB1->getWidth(), pB1->getHeight() + pB4->getHeight());

	int nRow, nColumn, nRover; // multipy B2 and B3 and add B1 writing results inplace
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = pB1->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow, nColumn) = pB1->valueAt(nRow, nColumn);
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB2->valueAt(nRow, nRover) * pB3->valueAt(nRover, nColumn);
		}

	// simple multiply of B4 and B5
	for (nRow = pB4->getHeight(); --nRow >= 0; )
		for (nColumn = pB5->getWidth(); --nColumn >= 0; )
			for (nRover = pB4->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow + pB1->getHeight(), nColumn) +=
					pB4->valueAt(nRow, nRover) * pB5->valueAt(nRover, nColumn);

	return pB;
}

CScatBlockLP* CScatMatrixLP::computeNormalBlock1x2F(const CScatBlock *pB1,
	const CScatBlock *pB2, const CScatBlock *pB3, const CScatBlock *pB4, const CScatBlock *pB5)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL && pB5 != NULL);
	// first expression is valid
	assert(pB1->getWidth() == pB2->getHeight());
	// second expression is valid
	assert(pB3->getWidth() == pB5->getWidth());
	assert(pB3->getHeight() == pB4->getHeight());
	assert(pB4->getWidth() == pB5->getHeight());
	// both expressions
	assert(pB2->getWidth() == pB3->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(false,
		pB3->getWidth(), pB1->getHeight() + pB3->getHeight());

	int nRow, nColumn, nRover; // simple multiply of B1 and B2
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = pB2->getWidth(); --nColumn >= 0; )
			for (nRover = pB1->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB1->valueAt(nRow, nRover) * pB2->valueAt(nRover, nColumn);

	// multipy B4 and B5 and add B3 writing results inplace
	for (nRow = pB3->getHeight(); --nRow >= 0; )
		for (nColumn = pB3->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow + pB1->getHeight(), nColumn) = pB3->valueAt(nRow, nColumn);
			for (nRover = pB4->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow + pB1->getHeight(), nColumn) +=
					pB4->valueAt(nRow, nRover) * pB5->valueAt(nRover, nColumn);
		}

	return pB;
}

// some as above - the formula is:
// | B1 + B2 * B3,	B2 * B4			|
// | B5 * B6,		B7 + B5 * B8	|
CScatBlockLP* CScatMatrixLP::computeNormalBlock2x2(
	const CScatBlock *pB1, const CScatBlock *pB2, const CScatBlock *pB3,
	const CScatBlock *pB4, const CScatBlock *pB5, const CScatBlock *pB6,
	const CScatBlock *pB7, const CScatBlock *pB8)
{
	assert(pB1 != NULL && pB2 != NULL && pB3 != NULL && pB4 != NULL);
	assert(pB5 != NULL && pB6 != NULL && pB7 != NULL && pB8 != NULL);
	// first expression is valid
	assert(pB1->getWidth() == pB3->getWidth());
	assert(pB1->getHeight() == pB2->getHeight());
	assert(pB2->getWidth() == pB3->getHeight());
	// second expression is valid
	assert(pB2->getWidth() == pB4->getHeight());
	// third expression
	assert(pB5->getWidth() == pB6->getHeight());
	// forth expression is valid
	assert(pB7->getWidth() == pB8->getWidth());
	assert(pB7->getHeight() == pB5->getHeight());
	assert(pB5->getWidth() == pB8->getHeight());
	// all expressions
	assert(pB1->getWidth() == pB6->getWidth());
	assert(pB7->getWidth() == pB4->getWidth());

	// allocate target block
	CScatBlockLP *pB = new CScatBlockLP(false,
		pB1->getWidth() + pB7->getWidth(), pB1->getHeight() + pB7->getHeight());

	int nRow, nColumn, nRover; // multipy B2 and B3 and add B1 writing results inplace
	for (nRow = pB1->getHeight(); --nRow >= 0; )
		for (nColumn = pB1->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow, nColumn) = pB1->valueAt(nRow, nColumn);
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn) +=
					pB2->valueAt(nRow, nRover) * pB3->valueAt(nRover, nColumn);
		}

	// simple multiply of B2 and B4
	for (nRow = pB2->getHeight(); --nRow >= 0; )
		for (nColumn = pB4->getWidth(); --nColumn >= 0; )
			for (nRover = pB2->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow, nColumn + pB1->getWidth()) +=
					pB2->valueAt(nRow, nRover) * pB4->valueAt(nRover, nColumn);

	// simple multiply of B5 and B6
	for (nRow = pB5->getHeight(); --nRow >= 0; )
		for (nColumn = pB6->getWidth(); --nColumn >= 0; )
			for (nRover = pB5->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow + pB1->getHeight(), nColumn) +=
					pB5->valueAt(nRow, nRover) * pB6->valueAt(nRover, nColumn);

	// multipy B5 and B8 and add B7 writing results inplace
	for (nRow = pB7->getHeight(); --nRow >= 0; )
		for (nColumn = pB7->getWidth(); --nColumn >= 0; )
		{
			pB->valueAt(nRow + pB1->getHeight(), nColumn + pB1->getWidth()) =
				pB7->valueAt(nRow, nColumn);
			for (nRover = pB5->getWidth(); --nRover >= 0; )
				pB->valueAt(nRow + pB1->getHeight(), nColumn + pB1->getWidth()) +=
					pB5->valueAt(nRow, nRover) * pB8->valueAt(nRover, nColumn);
		}

	return pB;
}

/****************************************************************************
* CScatMatrixLP - private class for LP matrix								*
****************************************************************************/

/****************************************************************************
* Overrides																	*
****************************************************************************/
int CScatMatrixLPT::getWidth() const
{
	// switch Transform and return width or height
	switch (m_nTransform)
	{
	case CTransform::eTransform_r180:
	case CTransform::eTransform_fr0:
	case CTransform::eTransform_fr180:
	return m_pOwner->getWidth();
	case CTransform::eTransform_r90:
	case CTransform::eTransform_r270:
	case CTransform::eTransform_fr90:
	case CTransform::eTransform_fr270:
	return m_pOwner->getHeight();
	default: // invalid
	assert(FALSE);
	}
	// invalid
	return -1;
}

int CScatMatrixLPT::getHeight() const
{
	// switch Transform and return width or height
	switch (m_nTransform)
	{
	case CTransform::eTransform_r180:
	case CTransform::eTransform_fr0:
	case CTransform::eTransform_fr180:
	return m_pOwner->getHeight();
	case CTransform::eTransform_r90:
	case CTransform::eTransform_r270:
	case CTransform::eTransform_fr90:
	case CTransform::eTransform_fr270:
	return m_pOwner->getWidth();
	default: // invalid
	assert(FALSE);
	}
	// invalid
	return -1;
}

const CScatMatrix* CScatMatrixLPT::getLOrT() const
{
	// TODO: NOT IMPLEMENTED YET
	assert(FALSE);
	return NULL;
}

const CScatMatrix* CScatMatrixLPT::getROrB() const
{
	// TODO: NOT IMPLEMENTED YET
	assert(FALSE);
	return NULL;
}

const CScatBlock* CScatMatrixLPT::getBlock(eBlock nBlock, bool bCreateIfNotThere) const
{
	// TODO: NOT IMPLEMENTED YET
	assert(FALSE);
	return NULL;
}

const CScatMatrix* CScatMatrixLPT::getTransformed(
	CTransform::eTransform nTransform, bool bCreateIfNotThere) const
{
	// validate params
	assert(nTransform >= CTransform::eTransform_r0 && nTransform < CTransform::eTransform_count);

	// now combine trasfroms to get the real one
	CTransform::eTransform nCombinedTransform =
		CTransform::combineTransforms(m_nTransform, nTransform);

	// return matrix reference from owner
	return m_pOwner->getTransformed(nCombinedTransform, bCreateIfNotThere);
}

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatMatrixLPT::CScatMatrixLPT(CScatMatrixLP *pOwner, CTransform::eTransform nTransform) :
	m_pOwner(pOwner),
	m_nTransform(nTransform)
{
	// should never be null transformed
	assert(m_nTransform != CTransform::eTransform_r0);
}

/****************************************************************************
* End of file																*
****************************************************************************/
